
local bw_scatter=500
local bw_lpoly=400
local max=2000
local nq=0
local rhs = "hours"

*************************************************************************************
*************************************************************************************
*************************************************************************************
// DISCONTINUITY PLOTS 
*************************************************************************************
*************************************************************************************
*************************************************************************************
use "$datapath/001_nlsy_child_mother_cog.dta", clear

** Delete top bunching
drop if hours == 2080

keep if hours<=2400
egen x=cut(`rhs'), at(0 1(200)2400)
egen tag=tag(x)
label variable x "Average hours working during three first years of child"
local xtitle "Average hours working during three first years of child"

label variable afqt_mom "Mother's AFQT Score"
label variable spouse_atbirth "Spouse is Present at Birth"

local list="afqt_mom spouse_atbirth"


foreach var of varlist `list' {
		egen M`var'=mean(`var'), by(x)
		egen S`var'=sd(`var'), by(x)	
		egen N`var'=count(`var'), by(x)
		gen M`var'_ub=M`var'+(1.96*S`var'/sqrt(N`var'))
		gen M`var'_lb=M`var'-(1.96*S`var'/sqrt(N`var'))	

		lpoly  `var' `rhs' if `rhs'!=0, generate(x`var' s`var') nograph noscatter se(`var'_se) degree(1) bwidth(`bw_lpoly') n(500)
		reg `var' `rhs' if `rhs'!=0
		predict s`var'_lin
		predict s`var'_lin_se, stdp
		
		gen upper_ci = s`var' + 1.96*`var'_se
		gen lower_ci = s`var' - 1.96* `var'_se
		
		* to keep same y axis scale as child's cognitive
		local yscale ""
		local ylabel "ylabel(`ylabel_`var'')"
		if "`var'"=="Res_cog" {
			local yscale "yscale(range(-0.2 0.4))"
			local ylabel "ylabel(-0.2(0.2)0.4)"
		}

		twoway rarea lower_ci upper_ci  x`var' if x`var'!=0, color(gs10%90 ) ///
				|| line s`var' x`var' if x`var'!=0, lcolor(black) ///
				|| scatter M`var' x if tag==1 & x==0, mcolor(black) msymbol(oh)  ///
				|| rcap M`var'_ub M`var'_lb x if tag==1 & x==0, lcolor(black) ///
				scheme(s1mono) legend(off)  xlabel(0(`bw_scatter')`max') ///
				ytitle(`: variable label `var'', size(large))  xlabel(0(500)2500) xtitle(`xtitle', size(large)) xscale(range(0 2500)) ///
				title("") note("") `yscale' `ylabel'
		graph export "$figpath/002_Disc_X`rhs'_Z`var'_`bw_lpoly'_onlyAtZero.pdf", replace	
	
		scalar drop _all
		drop upper_ci lower_ci M`var' S`var' N`var'
}
